Model of surface instabilities induced by stress 
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We propose a model based on a Ginzburg-Landau approach to study a strain relief mechanism 
at a free interface of a non-hydrostatically stressed solid, commonly observed in thin-film growth. 
The evolving instability, known as the Grinfeld instability, is studied numerically in two and three 
dimensions. Inherent in the description is the proper treatment of nonlinearities. We find these 
nonlinearities can lead to competitive coarsening of interfacial structures, corresponding to different 
wavenumbers, as strain is relieved. We suggest ways to experimentally measure this coarsening. 
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Elastic effects can strongly influence the morphology 
of materials, and hence influence material properties. If 
noncquilibrium elastic energies build up, there are dif- 
ferent ways for solids to release that energy. One is by 
plastic deformation, involving dislocations, another is by 
elastic deformation, which is commonly seen in thin-film 
growth. A non-hydrostatically strained solid which is in 
contact with its own melt or vapor can partially release 
its elastic energy by a morphological instability at the 
interface. This strain relief mechanism gives rise to what 
appears to be a buckling of the surface into trenches, or 
islands, of a particular spacing. It was first predicted 
by Asaro and Tiller 0|. Experimentally, it has been ob- 
served and studied by Torii and Balibar who strained 
He'* crystals non-hydrostatically as well by Berrebar et 
al. Q in polymer crystals. Furthermore, it is often 
associated with the dislocation-free Stranski-Krastanov 
growth mode (also called island-on-layer mode) of epi- 
taxially grown thin films as being observed for Ge/Si 
InGaAs/GaAs @ and InGaAs/InP Since the inde- 
pendent rediscovery of the instability by Grinfeld j^] and 
Srolovitz 1^, it is often referred to as the Grinfeld insta- 
bility. 

Several approaches have been employed to study the 
instability. They are either based on static energy min- 
imization calculations by a variational principle 0, or 
on a dynamical interface equation which describes mass 
transport, mainly surface diffusion, under the influence 
of the chemical potential which comprises surface free en- 
ergy and elastic energy [8[-[l4|. Linear stability analysis 
[p|, p^|JTl] ] predicts conditions for the onset of instability. 
Spencer and Meiron and Yang and Srolovitz Q 
studied the nonlinear evolution numerically, whereby the 
surface profile evolved to smooth flat peaks with sharp 
deep grooves. These studies have been limited to dimen- 
sion d = 2: Within the interface formulation, unphysical 
sharp cusps form within the grooves, leading to numer- 
ical stability problems |^3|. To make connection to ex- 
periments, one requires a model which comprises a full 
nonlinear description, and which can be used in three 



dimensions. 

In this paper, we present a Ginzburg-Landau phase- 
field model of the phenomena. An order parameter field 
(j){f) determines whether one is in a hard solid phase, 
which supports shear, or a soft disordered phase, here- 
after called the hquid phase, which does not. The posi- 
tion of the interface coincides with the rapid variation of 
this field. Such an approach has been applied successfully 
to other moving-boundary-value problems, such as phase 
segregation and crystal growth [|l^ . Indeed, our model is 
numerically robust, can be implemented in three dimen- 
sions, and is readily generalizable. We show below that 
we recover the Grinfeld instability in linear and highly 
nonlinear regimes. We furthermore probe the transient 
dynamics during the morphological instability, finding 
that competitive coarsening of interface structures takes 
place. We suggest ways to measure this experimentally. 

The physical mechanism for the stress-driven morpho- 
logical instability can be understood easily. A stressed 
solid can partially relieve its stress by differentially mov- 
ing material from valleys to hills, buckling at a particular 
wavenumber. In the less constrained peaks, lateral relax- 
ation occurs, unlike in the more constrained valleys. The 
resulting stress gradient drives the instability by creating 
deeper valleys, thereby increasing the stress gradient, and 
sustaining the growth of the perturbation. At sufficiently 
small length scales, capillarity prevents the formation of 
sharp cusps. 

The model we propose is based on a Ginzburg-Landau 
approach in which the elastic strain is a subsidiary tensor 
variable coupled to a nonconserved scalar order param- 
eter. This approach is related to that of Onuki Jl6| , p7t , 
Onuki and Nishimori , and Sagui, Somoza, and Desai 
]l9t , which was used to analyze elastic effects in phase- 
separating alloys The coarse-grained Ginzburg- 
Landau free energy is: 

.F(0,u.,)=^j/(0,?.,,) + ^|V0p], (1) 

where integration over r is indicated by the subscript on 
the integral, Uij = ^{dui/dxj + duj/dxi) is the strain 
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and Ui is the displacement field. The bulk free energy 
density f{(l),Uij) is given by: 



is Ai 



0, the solid will be stressed whereas the liq- 
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where the first term describes a three- well potential with 
= being the liquid and = ±1 the solid phase, en- 
suring that the liquid-solid phase transition is first order. 
The potential depths are determined by the model pa- 
rameter a, which together with the parameter / being 
proportional to the surface tension, determines the inter- 
facial thickness. The second term shifts the energy, so 
that, for constant elastic coefficients, solid and liquid are 
at coexistence. The convenient choice g{(j)) — ^cjP — jcf)^ 
guarantees ||l^ that both bulk phases keep their equilib- 
rium values = (hquid) and = ±1 (solid). The cou- 
pling constant e is related to the externally applied stress. 
The trace of the strain tensor is V • m, and /ei(0,Uy ) is 
the isotropic elastic free energy [pT|]: 
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where k is the compressibility, and /i is the shear modu- 
lus in the solid phase alone. By construction, the shear 
modulus in the soft liquid phase is zero, whereas it stays 
nonzero and constant in the hard solid phase. Since the 
the solid phase supports shear, whereas the liquid phase 
does not, our phase-field order parameter has a transpar- 
ent meaning in the context of the liquid-solid transition. 

It is reasonable to suppose that the elastic field relaxes 
much faster than (p. Then the elastic field can be solved in 
terms of the order parameter using the condition of local 
mechanical equilibrium: 5J- /5ui = j(Jij = 0, where a 
summation convention over repeated indices is implicit. 
The stress tensor, o-y = jbuij, is then given by 



(7y = (£5(0) -I- • w)% -f 2/ig(0)(My - -^V • u). (4) 

The solution of this to first order in the shear modulus 
is: 



uid is stress- free. For a flat surface = 0(?;), the solu- 



tion of Eq. (o) in two dimensions is Ux 



V • w = TrA g(r) 



(5) 



G(r,f')V: V' W')M,0\r")g{r'% 



where g{r) = g(0(r)). 



V,Mj = A,, ~ {e/K)^, V, / G(f, r')g{f'), (6) 

J r' 

V^G{r,r') = (5(r-r'), and My (r,r") = Vi Vj G(r,f')- 
{Sij /d)S{r — r '). In the absence of external strain, that 



and 



Uyy{y) = —{e/K)g{y). Therefore, the solid will be uni- 
axially strained with e determining the strength of that 
strain. 

The elastic field can now be expressed in terms of the 
order parameter. Substituting the solution for the strain 
field gives the free energy in terms of alone. The long- 
range character of the elastic field appears through M. 
Assuming relaxational dynamics, the equation of motion 
is given by: 
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/'(0) 



fV^c^ + ^l^g'{<l>) hi<j>) 



with r being the mobility and 



(7) 



M0)= 2 / / [G(f,r')V', V;-M,,(f',f")+ (8) 

+ M.,ir,r') M,,{r',r"Mr') g{r")- 

Rescaling length and time scales r ^ f/X where A is 
a characteristic length scale, such as the wavelength of 
the perturbation, t — > tT/X^, rescales the parameters to 
(3 ~ rX^/Ta, e — l^/a/X and c = [laeji^. We obtain a 
dimcnsionless equation of motion: 



dt 



-/?[/'(0)-e2v20 + c5'(0) 



(9) 



with three parameters /3, e, and c, giving the mobility, 
capillarity, and shear strength, respectively. 

Numerical simulations on a discrete lattice were per- 
formed in two and three dimensions. Euler's method was 
used for the integration in time. The Green function was 
solved in Fourier space. For all simulations presented 
here the mesh size Ax = 0.01 or 0.005, the time step 
At = 0.1 or 0.05, /? = 1.0, and e = 0.01. This choice of 
Ax and e guarantees that the surface is resolved by at 
least 8 points. The parameter set, (L^., Lj,, L^, YJ), c) will 
be specified below, where Iq gives the initial amplitude 
of the surface. Length scales will be measured in units of 
Ax. Periodic boundary conditions were employed in all 
direction. Thus, the solid was in contact with its liquid 
phase at the bottom and at the top. It was ensured that 
the solid was sufficiently thick that the interfaces at the 
top and bottom acted independently. 

A numerical linear stability analysis was performed 
in two dimensions. The system was prepared with a 
small amplitude sinusoidal surface profile Y[x,t = 0) = 
Yosin{qx), where q is wavenumber, and its subsequent 
evolution was monitored. We found that the growth of 
the amplitude of the Fourier modes was initially indepen- 
dent and exponential, obeying exTp{uj(q)t), followed by 
slower constant velocity growth. The fitted dispersion 
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ijj{q) is consistent with uo = Aq — Bq^, where ^ « 0.2 
and B « 25. See Fig. 1. Perturbations with wavenum- 
ber larger than a critical wavenumber are stabilized by 
surface tension, whereas wavenumbers smaller than the 
critical wavenumber are unstable, therefore being a long 
wavelength instability. The flat interface however is sta- 
ble. This agrees with the linear stability analysis car- 
ried out by Srolovitz jsj for the case where evaporation- 
condensation is the material transport mechanism, which 
is appropriate for our model. 
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FIG. 1. Early-time amplitude of growing mode uj plotted 
as ui/q versus wavenumber q to show that dispersion relation 
is consistent with uj — Aq — Bq^ . (L^ = Ly = 256, c = 5.6, 
Yq = 4, 1000 time steps). Lower inset: uj vs. q. Upper inset: 
Time evolution of a configurations in two dimensions sampled 
every 200 time steps. [Lx = Ly — 512, Yq — 12.1, q = 4/Lx, 
c — 5.6, 2000 time steps). Linear theory describes only the 
initial stages of the instability before asymmetry becomes ap- 
parent. 

Linear stability analysis predicts only the condition of 
onset of instability. To study the later-stage morphol- 
ogy, a complete nonlinear description has to be employed. 
One advantage of the phase-field description is that non- 
linearities are taken into account implicitly. A typical 
set of configurations is shown in Fig. 1. The nonlinear 
effect gives rise to a clear asymmetry between peaks and 
valleys, wherein deep grooves appear in the valleys. This 
behavior has been observed experimentally, as well as in 
previous theoretical studies P,p^-p^. Unlike previous 
studies, no numerical instabilities limit the study of the 
formation of the grooves here. It is interesting to note 
that in the early stages of the instability we can fit the 
interfacial profile with a simple fmiction K = 
where the curvature K = Y" (x) / {1 + Y' (x)'^)'^/^ is a low- 
order polynomial function of the height Y{x) of the in- 
terface. 

Experimentally, random fluctuations in the interface 
will give rise to the competitive growth of different struc- 
tures corresponding to different wavenumbers. To study 



this, we prepared the system with an interfacial profile 
consisting of a superposition of p linearly unstable modes, 
Y{x) ~ Yq X)r=i cos{qiX + (pi) with < and (p being 
a uniformly distributed random variable in the interval 
[0,2^]. 
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FIG. 2. Structure factor at equal time intervals. 

Solid lines show the structure factor derived from a linear 
Cahn-Hillard-type theory, which only describes the data for 
early times. For later times, structure factor is consistent with 
scaling, shown in inset. 

We did 100 runs over 500 time steps of a two- 
dimensional system with 100 unstable modes, where {L^, 
Ly, L^, Fo, c) = (1024, 512, 0, 0.24, 11). Figure 2 shows 
the Fourier transform of the equal-time height-height cor- 
relation function, which we shall call the structure factor 
S{q,t). For early times, there is a strong similarity be- 
tween this behavior and early-stage spinodal decomposi- 
tion in long-range force systems ; we show the results 
of a linear Cahn-Hilliard-type theory of the modes in the 
figure as well. Note that the structure factor vanishes for 
q ^ due to elasticity, not a conservation law. For later 
times, when the linear theory no longer describes the 
data, coarsening is evident: The location of the peak of 
the structure factor (7max(0 moves to smaller wavenum- 
bers, as the peak height increases and sharpens. The 
peak height follows S'(qmax,^) ~ t", where a ~ 2, while 
the peak width sharpens with time as w ~ t~'^ , where 
7 « 0.5. The former dependence is due to the total inter- 
face length increasing linearly with time for any unstable 
wavenumber. The latter dependence is due to competi- 
tive ordering between different wavenumbers, analogous 
to phase ordering. Within the accuracy of our study, 
we find that the structure factor shows scale invariance: 
S{q,t)/ S{qnia.x,t) — S*{q*), where the scaled wave num- 
ber q* ^{q- gmax)/?«- See Fig. 2. Fitting to S* ~ {q*Y 
and S* ^ (l/g*)'^, for small and large q* respectively, 
gives (5 ~ 1 — 2, and ^ 5 ~ 6. 

Although these results were obtained in two dimen- 
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sions, we expect qualitatively similar results in three di- 
mensions. To show this, we simulated a system with 
Lx = Ly = Lz = 128, with z being the direction normal 
to the surface. 
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FIG. 3. Typical configuration in three dimen- 

sions after 150 time steps, during the coarsening regime. 
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128, Yo = l,c= 16.7). 



Starting with a small amplitude sinusoidal perturba- 
tion in X, trenches with sharp deep grooves form, while 
a small amplitude sinusoidal perturbation in the x and 
y directions resulted in islands. The instability is quali- 
tatively the same as in two dimensions. Starting with a 
superposition of unstable modes, coarsening was again 
observed. Figure 3 shows the interfacial profile while 
coarsening is taking place. We expect that our results on 
transient coarsening phenomena can be observed through 
microscopy or by x-ray diffraction |^^ . 

To conclude, our model recovers the main features of 
the Grinfeld instability. Our description can be easily 
extended. Anisotropic effects can included through the 
surface tension, the elastic coefficients, or the external 
stress. The effect of phase separation or of impurities can 
be studied by coupling an additional field to the phase 
field. Instead of evaporation-condensation, surface dif- 
fusion can be chosen as the material transport mecha- 
nism, and, in addition, the influence of a constant flux 
can be studied. Finally, we note that in some cases the 
stress field at the groove tip can become so high that 
dislocations can be nucleated [^|p^. To study this, we 
are presently extending our model by coupling the phase 
field to a dislocation density field. 
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